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ABSTRACT 

Classical </> 4 theory in weak and strong thermal gradients is studied on the lattice in (1+1) 
dimensions. The steady state physics of the theory is investigated from first principles and 
classified into dynamical regimes. We derive the bulk properties associated with thermal trans- 
\Q . port, and explore in detail the non-equilibrium statistical mechanics of the theory as well as 

connections to equilibrium and irreversible thermodynamics. Linear response predictions are 
found to be valid for systems quite far from equilibrium and are seen to eventually break down 
simultaneously with local equilibrium. 
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I. INTRODUCTION 

The physics of non-equilibrium systems includes a broad class of phenomena, such as the physics of steady states, 
relaxation and dynamics far from equilibrium. Dynamical processes which range from those in the early universe, to 
ultra-relativistic heavy ion collisions and the formation of the quark-gluon plasma all involve non-equilibrium physics 
in an essential manner, perhaps requiring an understanding of the physics beyond linear response. Physics of steady 
■ states contains interesting non-equilibrium phenomena, such as transport, spatially varying observables, and is also 
important to the study of systems where the time scales for local equilibration are smaller than macroscopic time 
scales which might describe some global evolution of the system. Non-equilibrium steady states can be realized in 
many ways, such as the placement of a systems in a thermal gradient, or in an environment which provides shearing, 
pressure gradients and so forth. In this article, we will examine the statistical mechanics of the non-equilibrium steady 
states of a classical field theory in thermal gradients. This will allow us to understand the behavior of the theory 
under these non-equilibrium conditions and to consider problems related to the range or validity of local equilibrium, 
linear response, equilibrium thermodynamics and statistical mechanics. In particular, the thermal gradients we study 
can be quite strong and in such situations, it is natural to ask whether the system always relaxes to local equilibrium. 

When systems are near equilibrium, one expects linear response to provide a description of the transport coefficients. 
However, there is no means to address its regime of validity within the theory itself. Furthermore, even in this regime, 
in low dimensions, d < 2, it has been argued using kinetic theory that linear response often does not hold KJ; the 
Green-Kubo autocorrelation functions are expected to behave as t~ d / 2 , leading to divergent transport coefficients. 
Such divergences of the Green-Kubo integral have been observed in certain low dimensional systems such as the FPU 
model 101 or the diatomic Toda lattice 101 , and seem to be endemic to Hamiltonians which conserve total momentum. 
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Other studies have also found thermal transport in the linear regime to diverge |y,|5| or have focused on somewhat 
exotic models |]]7j , while strong thermal gradients have been studied in cellular automata || . 

In this work, we attempt to address many of the basic questions regarding the non-equilibrium properties of the 
(f> 4 theory on the lattice in (1+1) dimensions, by studying the model both near and far from equilibrium. We choose 
the <f> theory since it is a prototypical model which appears in a variety of contexts, including particle physics. From 
the outset, we should point out that our work has two limitations, namely that it is classical and that it is a lattice 
field theory. On the other hand, we make no further approximations and we analyze the model from first principles 
without any dynamical assumptions. This will allow us to answer interesting physical questions that cannot yet be 
addressed in the full quantum case. Furthermore, the approach we adopt can be generalized to other classical lattice 
field theories in a straightforward manner. Our main objective is to develop a comprehensive understanding of the 
underlying dynamics of the scalar field theory in thermal gradients and to lay the ground work for further analysis. 
As such, we shall provide the necessary details of our methods for further work and in such a manner so that can be 
easily generalized to other models. 

Classical field theory is relevant to high temperature behavior of quantum field theories: For instance, recently, it was 
used to derive properties of the standard electroweak model at finite temperature |j|,[l0| . Properties of the quantum </> 4 
theory in equilibrium has also been studied previously |ll|Jt2] ]. However, the relation between the physical quantities 
in the classical lattice field theory and the quantum theory is far from trivial and it is beyond the scope of this work. 
The dynamics of classical field theories is of interest on its own right, which in our case can be thought of as the 
dynamics of an anharmonic chain. While we discuss the results of the (1+1) dimensional theory in this work, it is 
not an essential limitation of our model or approach; we do find that the basic physics understanding developed in 
(1+1) dimensions carries over to results in (3+1) dimensions, although the latter will be discussed elsewhere. As such, 
our analysis in (1+1) dimensions is not specific to 1-d systems. It should be noted that classical 4> A theory has been 
considered in various contexts in the past: The Lyapunov spectra has been computed in the microcanonical ensemble 
in massless ]l3| and massive jl4| models. Equilibration of the model was studied in [ flEf . 

In addition to elucidating the physics underlying some of the non equilibrium phenomena of field theories, we 
believe that our results shed light on the nature of non-equilibrium statistical ensembles. The extension of the Gibbs 
ensemble to the non-equilibrium steady state remains unanswered. However, limited results to date suggest that the 
theory is far from trivial, including divergent Gibbs entropy and singular steady-state measures Hi. While some 
approaches exist, such as maximum entropy or projection techniques to construct non-equilibrium operators |l6| ], 
assumptions must be made about the non-equilibrium state in order to compute its properties. By using classical 
field theory as our starting point, we can use existing techniques to construct non-equilibrium steady states without 
assumptions on the dynamics of the model which are symptomatic to other approaches. A seemingly simple question 
concerns the thermal profile T{x) which develops in a system when it is in a thermal gradient. In other approaches, 
often some form is assumed for the profile T(x), which in principle should be dynamically obtained, as we shall 
do here. Furthermore, in our study, whether local equilibrium is achieved is not an assumption but is determined 
dynamically by the system. We will see that there are qualitative differences between the behavior of the system in 
the various regimes. These are characterized in Table 1. 

We would like to emphasize once again the motivation for this analysis: While the lattice model we work with 
does have well defined thermal transport, a single component scalar continuum field theory does not support thermal 
conductivity in the usual sense. It would be interesting to generalize to more complex theories with more than one 
conserved current, such as the two component scalar field theory, gauge theory with matter and so on. On the other 
hand, rather than introduce additional degrees of freedom and observables such as charge or matter density, we prefer 
first to focus on some important questions in the lattice theory. By focusing on the simplest lattice theory, we are 
able to clearly present an approach to the non-equilibrium statistical mechanics of classical lattice models from first 
principles, which can then be applied to other models like the ones mentioned above. Our results are also of interest 
to the statistical mechanics of many body systems in non-equilibrium. So as is, we do not attempt to make claims 
concerning the continuum limit but rather concentrate on elucidating the physics of the lattice theory near and far 
from equilibrium. 

In section 2, we describe the model we study and how we analyze the theory, particularly paying attention to 
the way the temperature boundary conditions are implemented. In section 3, we discuss thermal transport in our 
theory. In section 4, we analyze the equilibrium physics of the model and in section 5, the non-equilibrium physics. In 
particular, we study the thermal conductivity and its temperature dependence. We analyze the thermal profiles and 
establish that the linear response theory works even for visibly curved profiles, up to certain thermal gradients. We 
further examine the relations between the various physical quantities such as entropy, speed of sound, heat capacity 
and the thermal conductivity. We end with a discussion in section 6. 

Table 1: Behavior of the (jr theory under varying thermal gradients. 
Here, I is the mean free path as explained in section M. 
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Regime 



Properties 



Global Equilibrium (GE) 
Local Equilibrium I (LE-I) 



Local Equilibrium 1 (LE-E) 



Local Non-Equilibrium (LNE) 



VT = 0, f{v,4>)~exp[-H/T\. 

AT/T < 1; VT ^constant; 

Fourier's law holds globally; 

Agreement with linear response theory. 

IVT/T < 1/10; VT ^constant; 

Fourier's law holds locally; 

Small deviations from linear response theory; 

Existence of boundary temperature jumps. 

e\/T/T > 1/10; VT ^constant; 

Local equilibrium description inadequate; 

Definition of temperature ambiguous; 

No suitable definition for transport coefficients. 



II. THE MODEL 



We start with the (j) 4 Lagrangian (with the metric convention (—,+)), 



(1) 



We discretize and perform the rescaling 

4>x{t) — o.g<j>(x, t), t — t/a, x — x/a, m 2 — m 2 a 2 , (2) 
where a is the lattice spacing. We then obtain the corresponding Hamiltonian where the lattice spacing is scaled out 

1 



7 rf + (V0 l ) 2 +m 2 ^ + ^ 



(3) 



Here k = 1, 2, . . . L runs over all sites in the lattice, the lattice derivative is V(f>k = 4>k+i — 4>k- The resulting equations 
of motion are: 



<f>i = 7Tj, iii — (y 2 4>)i - m 2 ^ - (f>f. 
Here, we defined the lattice Laplacian as (V 2 ^)fc = cf>k+i — 2<t>k + <fik-i- 



(4) 



A. Finite Temperature Equilibrium Dynamics: VT = 

Starting from the microcanonical dynamics (^), we can develop a realization of the constant temperature dynamics 
using the global demons of Ref . . In this approach, auxiliary variables are added to the systems which dynamically 
emulate the presence of a heat bath. This type of dynamics, while not as rigorously understood, has been shown 
to converge much faster than the optimized hybrid Monte Carlo methods. Further, it does not suffer as much from 
critical slowing down |19|. 

When we are interested in studying the statistical properties of a system described by an action S(ip), where 
tp = ((pi, tf2, ■ ■ ■ ifin) are the degrees of freedom, we usually start with the definition of a statistical measure, such as 



fVn{<p) ~ exp[-S(<p)/T\V/i(<p). 



(5) 



Here, T>(i(ip) might include constraints in the dynamical space, as in the case of motion on curved manifolds, such as 
Lie groups. Because we know the measure, steady-state values of observables are readily determined. For an arbitrary 
observable O, we have 



(O) = - J T>p{<p) e- s ^' T 0, Z = j V^ip) 



(6) 
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While the approach we discuss is suited to general measures, in this article we use the measure T>/i(tp) = Dtp over 
the phase space, where tp will typically represent canonically conjugate coordinates and momenta, tp = (</>, 7r), and 
S(tp) is taken as a Hamiltonian, S(4>, tt) = H(4>, it). While the dynamics of the model may now be easily implemented 
using equations of motion, ([|) — often referred to as the molecular dynamics method — we would like to add finite 
temperature constraints to the equations. For this to happen, we must no longer evolve on the constant energy 
surface, so that S(tp) should no longer be conserved. The method we discuss here is reminiscent of Parisi and Wu's 
stochastic quantization [ p0[ , although the one adopted in our work is deterministic and time-reversal invariant. It 
is also a versatile approach in that it has been applied to systems with non-trivial measures, such as Lie algebraic 
Hamiltonians, equilibrium and non-equilibrium quantum systems, atomic clusters and molecules, magnetic materials 
and lattice models f2l]j . 

There are many formulations of this dynamics, initially motivated by the approaches of Nose and Hoover p2[ . 
Consider the following equations of motion for a thermostatted site labeled by k: 

dS dG(w k ) „. . dG'iw') . 

<f>k = TTfe, TTfe = -— F {n k ) - V F'(7T k ). (7) 

ocpk dwk dw k 

We have added two additional degrees of freedom, w k ,w k , which couple through forces indicated above. These extra 
degrees of freedom, so called 'demons', may be coupled either to the fields, <f) kl or to their 'momenta', n k . Here, 
we choose to couple them only to Tr k s in order to have the ability to apply thermostats locally at any one site. 
The microcanonical limit is recovered when these extra degrees of freedom are decoupled. In this extended space, 
tp = {4>i, iTi,w k ,w' k ) , we define a new action which is the old one plus additional terms for the demons: 



f(<j),ir,w,w') = exp 



Sfair) + J2 (G(w k ) + G'(w' k )) 



k : thermostatted 



IT . (8) 



In contrast to microcanonical dynamics, where the Hamiltonian is a constant of the motion, / is not preserved by 
the constant temperature dynamics. While the choice of forces in the equations of motion as well as that of / are 
seemingly arbitrary, steady state expectation values will be independent of these under reasonably general conditions. 
Consequently, they are chosen to optimize convergence of the physical variables. 

To find the dynamics associated with the demons u^u^, we simply require that / satisfy a continuity (Liouville) 
equation in the configuration space ip = (0.;, 7r.;, w k ,w' k ): 

0= »/ + V-W) (9) 

dt + ^ dtpi ■ w 

This is equivalent to requiring that the master equation, enforcing conservation of probability under evolution of the 
ensemble, be satisfied. By substituting the equations of motion into the continuity equation, and using the definition 
of /, we can derive solutions for w k ,w' k : 

771/ X rr.dF^k) ■ dF'(-K k ) 
w k = TT k F{ir k ) - T— , w' k = K k F (7T fe ) - T — . (10) 

dTT k dTT k 

By construction, this dynamics preserves the measure Eq. (||), so that time averages of observables on a given trajectory 
will converge to the configuration space average over the canonical measure. There is clearly some freedom in defining 
the dynamics; namely, the functions G(w),G'(w') and F(ir), F'(tt). The only restriction on G(w),G'(w') is that the 
measure Eq. (^) leads to a finite integral; in general the auxiliary variables w can have any desired measure. In 
practice, highly non-linear functions are impractical since they will require small integration time steps. For these 
reasons, it is convenient to take G(w), G'{w) to be fiw 2 /2 or ^'u> 4 /4, where fi, pi! are positive constants. The constant 
couplings fi, n' of the demons to the physical degree of freedom are in principle arbitrary as long as the phase space 
integration is finite. Choosing these couplings to be too weak will make the convergence slow or choosing them to be 
too strong will lead to small time steps in the evolution, so that the couplings are chosen to optimize the convergence 
of physical observables. A necessary condition for F(n), F' (n), on the other hand, is for it to be at least linear in 
its argument, the minimal requirement for the existence of the fluctuations in the phase space volume which allows 
for the exploration of the canonical measure. The precise relation to the fluctuations in a phase space volume V, or 
equivalently, the instantaneous phase space compressibility, can be found using the divergence theorem 

IdV _ f ^ ^ / dG{w k ) dF(ir k ) ( dG'{w' k ) dF'{-K k ) 



i Vip y { ^-^-^ + — v fc/ v ^ ) , (ii) 

V dt J v t-r 1 \ dw k d-K k dw' k dir k 
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where the index k runs over the thermostatted sites. In this paper, we do not explore the effect of different choices of 
G(w), G'(w), F(tt), F'(tt). Such studies have been done on various other systems 

Finally, we note that the linearized equations of motion are evolved by the stability matrix, d(pi/d(pj. The eigen- 
values of this matrix are the Lyapunov exponents for the system. Hence we have the relation that 

i i 

For the canonical ensemble, the Liouville equation gives ^2 = 0? while in steady state non-equilibrium systems, 

E a. • i). 

One specific realization of the finite temperature equilibrium dynamics we use couples the thermostats only at the 
endpoints of the systems, k = l,L. With the choice of G(w) = w 4 /4, G'(w') = w' 2 /2, F(tt) = vr/T, F'(tt) = n 3 /T, we 
obtain 



9k 




TTfc 


ill 







k = 1,2,... L 

(V 2 0) fc - m 2 fc - 03 fc = 2,3,...,L-l 

(V 2 0) fe - m 2 (bk -<t>%- w 3 k iTk/T - w' k ir 3 k /T k = 1, L 
/T-l, w' k =4/T-3n 2 k k=l,L. 



(13) 



Just thermostatting the two boundary points is sufficient to thermalize the entire system. We then can examine 
the interior system far from the boundaries to study the finite temperature theory. Either free or fixed boundary 
conditions were used for the field with no significant effect on the physics behavior of the theory. To compute 
observables, we use the fact that in this type of dynamics, the time averages converge to the ensemble average using 
the desired ensemble (||): 

O = Urn - (dt O(0(t), n(t)) = (0) EQ = . (14) 



B. Non-equilibrium Boundary Conditions: VT 7^ 

One of the problems immediately encountered in the study of non-equilibrium systems is the nature of the steady 
state statistical distribution. While equilibrium statistical mechanics is well understood, once we apply thermal 
gradients to the system, much less is known about the system. To set up the non-equilibrium molecular dynamics, 
we use the demons to thermostat the end-points of our system in the same way we did in the equilibrium simulation. 
The only difference is that we now control the two endpoint temperatures separately. A consequence will be that the 
phase space distribution, f{4>, tt, t), will evolve to a non-smooth function that describes the non-equilibrium steady 
state. 

We take the equations of motion at finite T and now introduce two temperatures T® and T® . The superscript is 
needed to distinguish the thermostatted temperatures from those measured just inside the system, which can suffer 
from boundary jumps which we will analyze in detail |23j |. One set of equations, similar to the equilibrium case in 
Eq. ([l^), we have used are 

4>k = TT/c fe = l,2, ...L 

f (V 2 0) fe - m 2 fe -(A k = 2,3, .... L - I 

^ = { (V 2 0)! - m 2 0! - 0? - wfm/T? - w'^l/n , , 

I (v 2 0)l - - <t>i - wi-K L m - «m [ ' 

m = Tif/T? - 1, w[ = 4/T? - 3tt 2 
w L = nl/T°-l, w< L = ni/T°-3irl. 

One can see that the thermostats are applied to the endpoints k — 1 and k = L. It should be noted that inside the 
boundaries, the dynamics of the system is that of the 4> 4 theory itself with no other degrees of freedom. We have also 
considered both variations of these thermostats, such as different forms of the interactions or increasing the number 
of sites at each end where we apply the demons. We will comment when these distinctions are relevant. 

The equations of motion are solved on a spatial grid, using two methods: fifth and sixth order Runge-Kutta, and 
leap-frog algorithms [Q. We used from 10 6 to 10 9 time steps of dt from 0.1 to 0.001, with observables being sampled 
every At = 20 ~ 100 dt. The lattice size was varied from L = 20 to 8000. 
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A consequence of the non-equilibrium steady state is that the measure becomes singular with respect to the Liouville 
measure |17|] . We do not rely explicitly on the singular nature of the non-equilibrium measure. Rather, we use it 
to interpret certain observables in the non-equilibrium state. To see this we start with the two main equations for 
f(<fi,ir,t); the continuity equation and the expression for the total derivative of a phase space valued function. If we 
denote the vector p to include all degrees of freedom, p> = (0, 7r, u>), we have 

df df . df 

By combining this equation with the continuity equation which holds both in equilibrium and in non-equilibrium, 
we derive 

d£ = _ydp 1 



dt dp 

Solving for /, we obtain 



fit) = HO) exp - 1 it £ _ /(„) exp - t £ (g) _ /(„) exp U £ A, . (18) 



In these steps we have replaced the time average with the non-equilibrium ensemble average. The average of the 
divergence of the equations of motion is nothing more than the sum of Lyapunov exponents. In non-equilibrium 
steady states, whether from thermal gradients, or shearing, and so forth, the sum of the exponents is observed to 
become negative, signaling the presence of a fractal dimension. In the steady state, 

/(*)t=£,°°- ( 19 ) 

As the continuity equation is satisfied, the allowed phase space volume must be shrinking onto a set of measure 
zero, with respect to the original measure. A consequence of the divergence of the distribution function in the 
non-equilibrium steady state is that Gibbs entropy will also diverge, 

S G = -(log/) - -oo. (20) 

Although we do not know how to properly define the fractal measure for our non-equilibrium steady state, we are 
able to compute non-equilibrium expectation values with respect to it using time averages: 

1 '* 



(0) NE = 0=\im- dtO(4>{t),-K{t)). (21) 

t-,oo t J 

We will verify this in the linear response regime where near-equilibrium results can be compared to thermal equilibrium 
predictions obtained using the linear response theory. 

Simulations of many-body systems in non-equilibrium steady states have found that the steady state measure is 
typically singular with respect to the original equilibrium measure. As one moves further from thermal equilibrium, 
the available phase space contracts onto an ergodic fractal: the accessible points are dense in the phase space, but 
fractal in nature. The resulting loss of dimension is related to the transport coefficient. One can see that in our 
dynamics as well since / satisfies the continuity equation, so that total probability is conserved, yet Eq. (19) is 
satisfied. This means that the accessible phase space volume is contracted onto a set of measure zero. This type of 
'dimensional loss' in the steady state can be more rigorously seen in low dimensional systems like the Lorentz gas. 

From the point of view of dynamical systems theory, it is possible to understand the steady state measures under 
certain special conditions, namely that the dynamics is hyperbolic or Anosov. Unfortunately the conditions for 
hyperbolicity are not satisfied for our system since the number of positive Lyapunov exponents can vary along a 
trajectory. Nevertheless, it is useful to note that for hyperbolic systems described by some flow, x(t) — StX, for time 
evolution operator St and initial point x, the Sinai-Ruelle-Bowcn theorem provides the existence of the steady state 
measure, denoted /j. S rb- It can be shown that for a continuous function f(x), there exists a unique measure fJ, SRB 
such that 

^\fdt>f( X {t))= ! f lf^ B . (22) 
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This measure is known to be fractal for various non-equilibrium systems, and can be explicitly constructed for certain 
maps such as the modified Baker's map JIt]]. Here a basis of fractal Takagi functions has been implemented. 

For the Lorentz gas, it was shown numerically under general conditions, and rigorously proven in the linear response 
regime, that the steady state measure is singular |P5| . Here one can explicitly see that the dimensional loss AD is 
proportional to the transport coefficient, which in this case is the electrical conductivity. While there is still some 
contention concerning the existence of singular measures in the extension of the Gibbs ensemble to the far from 
equilibrium steady state, the compelling evidence suggests that whatever its nature is, it is far from trivial. 



III. TRANSPORT 



One of the relevant observables is the stress-energy tensor, 



dC 

From the continuity equation = d^T' 11 ', we have 

d 3T 00 dT 0i 

° = d^ T " = dt + ~dx T ■ ^ 
We can identify the heat flux as T 0i . On the lattice and in one spatial dimension, 

T 01 (x) = -7r(aj)V^(x) -> T k 01 = -v k {<j> k+l - fe ). (25) 
Defining the lattice energy density consistently as 

T°°(x) = Itt 2 + i(V0) 2 + V(d>) - T°° = \*\ + i(0 fe+1 - fe ) 2 + V(cb k ), (26) 

we find that the discrete version of the continuity equation is satisfied: 

|^° + (r, 01 -T fe °i 1 )=0 . (27) 

This establishes that T k 01 is a constant in space and time in steady state. On the other hand, there is no way to 
satisfy the spatial component of the continuity equation, dtT k 01 + VT 11 = 0. The reason for this can be understood 
as follows; while the translation invariancc in the time direction is preserved, the translation invariance in the spatial 
direction has been lost due to the lattice discretization. 

It is interesting to contrast this with the so called FPU j3 model, which has divergent thermal conductivity in one 
spatial dimension. In this case the Hamiltonian is 



IE 

so that the heat flux reads 



Hfpu 

k 



pI + (ik+i - qkf + 7}(qk+i - q k ) 4 



(28) 



T 01 = -K k {q k+1 - q k ) (1 + P(q k+1 - q k f) . (29) 

If we think of this as originating from a continuum field theory, we would have the following (non-relativistic) La- 
grangian and heat flux: 

-4(t) 2 +KSM(S) 4 . ^-m(M%) 2 ) ■ 

It is straightforward to check that the heat current satisfies the continuity equation (pi)). In this case, the thermal 
conductivity diverges as L 2 / b ||. 

In the simulations we will define the temperature locally using the concept of an ideal gas thermometer, which 
measures the second moment of the momentum distribution. Specifically, at a site k (x — ka) we define in both 
equilibrium and non-equilibrium 

T(x)=T k = (4). (31) 

We will see that this definition is sensible in these regimes in the context of both statistical mechanics and thermo- 
dynamics. 
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A. Green-Kubo Approach 



The standard approach to near equilibrium or local equilibrium transport is to apply linear response theory or 
Green-Kubo formulas. In this approach, external fields are added to a Hamiltonian which generate the desired 
transport processes and expressions for the transport coefficients can be derived in the weak field limit. For thermal 
conduction, the arguments are somewhat heuristic since one does not have well defined external fields that produce 
heat flow. Nevertheless, the procedure empirically results in sensible expressions. The Green-Kubo approach expresses 
transport coefficients in terms of equilibrium autocorrelation functions. For the thermal conductivity, one has 

k (t) = ^>j™ dt J dx (r 01 (x,t)T 01 (x ,0)) EQ 

dt J2 (T 01 (x k: t)T 01 (x k ,,0)) EQ . (32) 

k,k' = l 

N(< L) is the number of sites in the region inside the boundaries used in the computation. We typically choose the 
region to be as large as possible while excluding the boundary effects. While this expression is expected to hold near 
equilibrium, the region of its validity cannot be determined within the linear response theory itself. 

Autocorrelation functions, such as Eq. (|32]), have been ar gue d to decay algebraically, rather than exponentially. 
Originally observed in early molecular dynamics simulations |26| , the velocity- velocity autocorrelation function was 
found to decay as (v(i)v(0)) EQ ~ t~ d l 2 up to times on the order of a few ten times the mean free time. This behavior, if 
it persists, leads to divergence of the diffusion constants for d < 3. Using kinetic theory and some general assumptions, 
it was later argued that this is a generic feature of dynamical systems for long times. These long-time power law tails 
cause divergences in transport coefficients in a large class of low dimensional systems including the FPU (3 model 
( p8| ) — (p0[). It is believed, however, that for theories without strict momentum conservation, such a divergence can 
be absent. This is the case for our lattice model due to our 'on-site' nature of the potential, as we shall see below. 




B. Non-Equilibrium Approach 

We also compute the thermal conductivity directly by constructing ensembles near equilibrium which lead to con- 
stant gradient thermal profiles. This is arguably more fool-proof than the Green-Kubo formula since no assumptions 
are necessary for the computation. Here we first confirm Fourier's law and then use it to compute the conductivity 
using 

<T) = (33) 
Vi (X) 

While Fourier's law is believed to be valid near equilibrium, what constitutes the 'linear regime' is not known. We 
will see that it eventually breaks down as we move far from equilibrium, while still in the steady state. It is not clear 
at this point how to make a sensible definition of the thermal conductivity. One approach would be to try to the 
characterize the 'nonlinear' response by incorporating the dependence of the conductivity on higher order terms in 
the derivatives, such as (VT) 3 , (V 2 T)(VT) and so on. We will try to clarify these issues below. 



IV. EQUILIBRIUM ENSEMBLE 

When the boundary temperatures are equal, we recover equilibrium physics for our system. In the remainder of 
the paper we specialize to the m 2 = case. 



A. Thermalization 



A test of the thermalization of the system is readily performed by studying the distribution functions of various 
quantities. If we take a single trajectory, (<fr(t),n(t)), and histogram 7Tfe(i) at t = 0, At, 2At,..., where At is some time 
step which ensures the points are reasonably well decorrelated, we will converge to the thermal distribution function 

/(7r fc ) ~ exp[-7r£/2T]. (34) 
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In Fig. [T] (left column), we show computed equilibrium distributions for the momenta and heat flux at the center of a 
lattice with L = 163. Here we have taken the endpoint temperatures to be = T§ = 1. The measured temperature 
in the middle is T = 0.995(8), and one can see that the measured momentum distribution (histogram, top left) agrees 
with the predicted value (solid, top left). In the bottom left we show the measured thermal distribution of T 01 , for 
which we do not have a theoretical prediction. Since we are in equilibrium, there is no heat flow, so (T 01 ) = 0, 
which is evident from the symmetric nature of the thermal distribution /(T 01 ). We have verified that when Tf = T§ , 
these boundary conditions reproduce the equilibrium canonical measure f e q(Tt,<p) cx exp[— H(jr, <p)/Ti] at all points, 
including the thermostatted sites. The measured thermal profiles satisfy T(x) — T®(= T$ ) for any x within numerical 
error over the temperature range we investigated: T = 0.01 to T = 10. Hence the boundary conditions ( |Te| ) do 
produce the desired physics for this theory. 

In our presentation of equilibrium and non-equilibrium steady state expectation values, we verify that stationary 
results are obtained. This is done by examining the time evolution of observables. For instance, if we want to measure 
the heat flow through the system, we measure at all sites and verify that the time averages converge to the same 
value. A typical result is shown in Fig. ^, where the average (T 01 ) is shown as a function of time at three sites: L/4, 
L/2 and 3L/4 for L — 800. One can see that the values eventually converge to the ensemble average. In this case the 
endpoint temperatures are distinct, so that there is a net heat flux. Equilibrium figures are similar, with convergence 
to (T 01 ) = 0. 

B. Auto-Correlation Functions 

We compute auto-correlation functions for the components of the stress-energy tensor, which includes the Green- 
Kubo formula for the thermal conductivity. We define the normalized correlation functions 

= W^W) [< r ^ T ^°)> - (^(O)) 2 ] • (35) 

In Fig. ||, we show the time dependence of | C(t) | for T 00 (dashes), T 01 (solid) and T 11 (dots) as a function of 
time, where the time is normalized by the mean free time r. (As we discuss below, r will be of order of the thermal 
conductivity.) One can see that these functions decay to half their initial values on the order of the mean free time. 
For larger times t > 10t, these functions oscillate about zero. This behavior is seen at all temperatures studied. 
To compute the thermal conductivity, we are interested in the auto-correlation function for T 01 : 

1 t* 

^(T,t) = -^ dtJ^i^i^t^ixk^O))^, K (T)=lim /c(T,t). (36) 

Jo k,k' 

By studying the convergence k(T, t) to the thermal conductivity, we can explore the problems associated with the 
long-time tails. In Fig. || we plot n(T,t) as a function of time, where t is normalized by the mean free time r. 
According to the predictions based on kinetic theory, in d = 1 we would expect k(T, t) ~ i 1 / 2 , leading to an infinite 
conductivity. At the temperatures shown (T = 1/50, 1/10, 1, 2), n(T,t) does apparently display behavior similar 
to i 1 / 2 (dashes) on time scales up to t ~ 10t, but on longer time scales, the results converge to well defined values. 
Consequently, divergences such as those found in the FPU model are not present here, and long-time tails are at best 
a transient aspect of this dynamics. 

The behavior of k(T) found from the Green-Kubo approach is summarized in Fig. || (crosses). This analysis will 
be continued when we discuss the direct measurements below. 



C. Speed of Sound 

To better understand the kinetic theory aspects of this finite temperature theory, we would like to know the thermal 
behavior of the 'speed of sound'. This can be estimated in a number of ways. A convenient approach is to use the 7~ 01 
auto-correlation function. The sound speed, c s , is defined here as how fast excitations travel through the system and 
is the relevant velocity for the transport theory in our particular model. We note, however, that strictly speaking, 
this is not 'sound' in the hydrodynamic sense. We define 

G(x,t;x ,to) = {T 01 (x ,t )T 01 (x,t)) EQ . (37) 

Instead of the time dependence, we consider the spatial dependence of G{x, t; Xq, to). Not only will it decorrelate in 
time, as we saw previously in the Green-Kubo integrals, but it will also decorrelate over space. A typical behavior of 
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G(x, t; Xo,to) is shown in Fig. |[ Here we choose xq — to be the center of the lattice, and to = 0. The temperature is 
T = 1/10 and L = 160. The autocorrelation function is shown for times t = 0, 30, 60, 90 where one can see the regions 
of high correlation separate. By measuring the rate at which the peaks separate, we can obtain an estimate for c s 
at that temperature. In Fig. 0, we plot the sound speed extracted in this manner against the temperature. There is 
little temperature dependence until T ~ 1/10, after which the speed begins to decrease with T, as one might naively 
expect. So although it shows some temperature dependence, c s is generally of order unity over this temperature range. 



V. NON-EQUILIBRIUM STEADY STATES 



By controlling the boundary temperatures, we can begin to develop an understanding of the physics of hot scalar 
field theory, as it moves increasingly further from thermal equilibrium |27| . We discuss the regimes we have found as 
we move away from equilibrium. 



A. Near Equilibrium: Ti ~ T 2 



In the near equilibrium limit, the physics is consistent with linear response and Fourier's law. Here (T 01 ) will 
develop a non-zero expectation value. As the temperatures T\ and T 2 begin to separate, a constant VT develops. 
This is illustrated in Fig. || for a lattice of L = 8000 where a small gradient is developed around a temperature 
T = 1/4 using boundary temperatures T® = 0.15 and T° = 0.35 (solid). A linear profile, expected from Fourier's 
law, is superimposed (dashes) . While such linear profiles in principle provide a direct measure of the conductivity by 
dividing the measured heat flux by the average gradient, we perform a more careful analysis by taking increasingly 
smaller gradients about the same average temperature, thereby explicitly verifying Fourier's law. This is shown in 
Fig. H for the temperature T = 1/4. All these points correspond to thermal profiles which are linear. As the boundary 
temperatures approach T = 1/4, we see that the heat flux also decreases. The slope, k(T) — — (T 01 )/VT, corresponds 
to the thermal conductivity at this temperature. 

The measurements of k(T) in the range 1/100 < T < 20 are summarized in Fig. |^ for both direct measurements 
(□) and Green-Kubo (x). We find that both methods agree and that the general behavior is well described by a 
power law: 

«(T) = A, 7 = 1-35(2), A = 2.83(4). (38) 



This type of behavior is reminiscent of lattice phonons at high temperature 28 1 . The fact that the direct measurements 
agree with the Green-Kubo integrals of the auto-correlation functions is sufficient to dispel the notion of asymptotic 
long-time tails in this system. 

Another important aspect to verify is that we have achieved a bulk limit in the values of k(T) reported. In Fig. 
[Toj , we show the dependence of the direct measurements of K on the size L of the system for several temperatures. 
The dashed lines are the predictions from the power law fit. One can see that a bulk limit is realized for relatively 
small lattices. On closer inspection, we find that the bulk limit is reached for smaller lattices when the temperature is 
higher. This can be understood as follows; the mean free path of the system is of the order of the thermal conductivity, 
as we shall see when we analyze the kinetic theory aspects of this problem. The bulk limit is reached for lattices 
larger than the mean free path, roughly speaking. 



B. Far from Equilibrium: Ti < T 2 : 



We now consider the behavior of the system as we move it further from equilibrium. One of the first characteristics 
to emerge is the development of curvature in the temperature profile. In Fig. hi] we plot a succession of steady 
state temperature profiles as we change the boundary temperatures from {T^T^) — (0.3,0.7) (dots) to (0.2,0.8) 
(dashes) and finally to (0.1,2) (solid). As the system moves away from the constant gradient profiles, it starts to 
feel the temperature dependence of the thermal conductivity, which decreases with increasing temperatures. As a 
consequence, the hot end of the system cannot conduct heat as well as to lower temperature end. (Of course the 
converse of this argument would hold if the power law behavior has 7 < 0, in which case the curvature would have 
an opposite sign.) 

Fig. |^, which was already mentioned when we discussed the convergence properties, displays the time evolution of 
the heat flux at three sites within a system thermostatted at T§) = (0.3, 0.7) (dots in Fig. |ll|). Regardless of how 
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far the system is from equilibrium, the heat flux must be independent of position in steady state, since there are no 
sources or sinks for heat inside the boundaries. We see that the values at the three sites converges to the same value 
in the steady state, namely, a constant non-equilibrium heat flux independent of x. 

Another property of the non-equilibrium system is the appearance of jumps in the temperature at the boundaries. 
This is illustrated in Fig. |l^; in Fig. [l^ (top) , we show a typical non-equilibrium thermal profile with curvature (solid 
— the dashes will be discussed below). In the lower panels we examine the low and high temperature ends. T® and 
T® are the temperatures enforced by our boundary temperatures. One can see that there is a difference between 
these temperatures and what one obtains by smoothly extrapolating the temperature profile near the edge. These are 
physical phenomena associated with the dynamics at the interface, which can be readily understood quantitatively, 
as discussed in section |y] D below. We will define the points obtained through smooth extrapolation by (T\,T2). If 
we now focus on the thermal profile away from the edges, we can understand the curvature in the temperature profile 
using Fourier's law. At or near equilibrium, we have found the power law behavior k(T) — AT 1 , as in Eq. (p8|). If 
we integrate Fourier's law and re-express the result in terms of extrapolated temperatures (Ti,T2), we find 



T(x) 



1-1- 



1-7 



7^1 
7 = 1 



(39) 



Some previous efforts to model the temperature profiles in systems exist [ p9| , although these were generally model 
fits and furthermore, a full understanding requires a description of the boundary effects, which up to now have been 
lacking. We have found that this formula provides an excellent description, but ultimately breaks down as we see 
below. In Fig. [l2|, the dashed line is the description given by Eq. ( |39| ) with 7 — 1.35, which offers a fit to within a 
few percent. In should be noted that 7 used here was obtained from systems at or near equilibrium as in Eq. d3§|), 
independently of the systems far from equilibrium. One can see that at the high temperature end there is a tendency 
to overshoot the measured behavior as in Fig. [l^ (bottom right). This analytic expression provides a sufficiently good 
description of the physics as we move away from equilibrium, until Fourier's law fails to hold locally. This includes the 
first three regimes listed in Table 1. Eventually, Fourier's law will break down locally as well, as we develop steady 
states which are locally non-equilibrium (LNE). 

The analytic understanding of the thermal profile allows us to understand the behavior of the heat flux for systems 
not too far from equilibrium. Using Fourier's law and Eq. ((3^), we have for 



(T 



01 \ 



.4 



NE 



L(l- 7 ) 



(Tl- 



T 



1-7 



) (7^1). 



(40) 



For 7=1, (T 01 )ne = (A/L) log(Xi/T2). When compared to measured heat flux using the extrapolated temperatures 
from the thermal profile, this formula is found to provide a very good description of the heat flux, typically to within 
a few percent. In the near equilibrium regime where the temperature profile is visibly linear, the boundary jumps 
vanish, and we expect it to behave as 



•01\0 



NE 



T 2 -T x 



(41) 



where the superscript denotes the linear response, or constant gradient limit. In general, these are related by 



)NE - (T 01 )°NE 



7(7- 



24 



11 f AT 



o 



AT 



T„ 



(42) 



where T av = (T 2 + T 1 )/2 and AT 



(T 2 — Ti). Note that these temperatures are the extrapolated ones and not 
In this way we see that the temperature profile will no longer be linear when 



the boundary temperatures, (T^T®) 
AT/T av ~ LVT(x)/T{x) ~ 1. 

In this regime of curved thermal profiles, we see some differences between the thermal distributions measured in 
equilibrium. In Fig. |l| we display on the right side the measured steady-state non-equilibrium statistical distributions 
for momenta (histogram, upper right) compared to the equilibrium thermal distribution expected at that temperature 
based on the notion of local equilibrium (solid). We find that as we more further from equilibrium, the momentum 
distributions are typically sharper than the gaussian expected based on the ideal gas thermometer. In the lower right, 
we show the non-equilibrium distribution of heat flux which develops the asymmetry needed to give a non-vanishing 
expectation value (T 01 ) < for T x ° < T 2 °. 
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An illustration of a system very far from equilibrium is the non-equilibrium steady state shown in Fig. [T^. Here the 
ratio of the endpoint temperatures is 100, and the thermal profile is shown on a log scale (solid). The fit of Eq. (|39| ) 
is shown by the dashed line. While we still measure temperature locally through the second moment, the measured 
distributions are no longer gaussian. In Fig. |lj (top) we show f(n) in the center of the system. In the bottom panel, 
the ratio of Eq. ( |39| ) to the measured behavior in Fig. [l3] is shown. Deviations up to 40% are evident on the low 
temperature end. Here, there is no accepted formalism to describe the dynamics and we are in the LNE regime of 
Table 1. We will return to this issue when we discuss local equilibrium in more detail in section M F. 



C. Entropy and Kinetic Theory 

While the Gibbs entropy for the theory diverges as it contracts onto a set of measure zero, not all measures of 
entropy behave in this manner. Irreversible thermodynamics provides a description of systems on scales much larger 
than the mean free path |3G] . When one has heat flow, one defines the local rate of entropy production as 



a(x) = (T t)1 ) NE 



1 



dx T(x) 



(43) 



i 2 (l-7W 



1-1^ 



1-7' 



2 r 




> 



Integrating this formula, we can compute the net rate of entropy production to be 

f L 

Sirr = dx a(x) 



(T 01 ) NE 



A 



£(7-1)1? 



1 1 

T 2 ~ Yi 



T 2 



(44) 



7-1' 



> 0. 



This expression can be simply interpreted as saying that the global entropy production rate is due to the difference in 
entropy production rates at the boundaries coming from the demons. On the other hand, irreversible thermodynamics 
predicts a constant local rate of entropy production, cr(x), which is at odds with the coarse-grained local Boltzmann 
entropy computed below. The latter calculation is microscopically based and does not rely on any hydrodynamic 
limit of the theory. 

If we envision the expansion of hot systems similar to RHIC collisions, the local entropy is an important quantity. 
From the statistical mechanics point of view, when one discusses the notion of entropy in a local frame, it is not Gibbs 
entropy of the entire systems that is needed. Hence we would like to consider Boltzmann's entropy. To understand the 
behavior of this entropy in increasingly non-equilibrium environments, we consider the n-body Boltzmann entropy, 

(n) 

defined through the n-body distribution functions f B . The best we can do, of course, is the coarse-grained limit of 
these quantities. In the 'local frame' at x — x' , we define them as integrals of the full phase space distribution / over 
all quantities except the arguments of f B . 

f < i\4>{x')Mx')) = f < i\ ( t> k ^ k ) 

f^i^x'l^x"),^'),^")) = f^\ct> k Ak+U^k + x) (45) 



The distributions f B are obtained by histogramming the corresponding degrees of freedom ( (0fe,7Tfe) 

(<j>k, 7Tfc, 4>k+ii TTfc+i) for f B ,...), and is readily constructed in equilibrium and non-equilibrium conditions. In the 
equilibrium or non-equilibrium steady states, we then compute 



for f<£\ 



o(l) _ 
"23 — 



s 



(2) 



d-K{ X )d4>{x) f B X) bg/W 
d-K{x)d<t>(x)d-K{x')d<t>( X ') f B 2) log/f } 



(46) 
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These satisfy the inequalities 

S G <---<^SW<LS£\ (47) 
To compute Sb we must coarse-grain the phase space. Hence we are actually computing the coarse-grained 1- and 

(k) 

2-body densities, which we denote . These are related by 

S B 1] ~ - J2 /a log /a - log( An k Afa ) , E /a = 1 ■ (48) 

We have computed these entropies and find that Sg^ does not shift noticeably from its equilibrium value regardless 
of how far the system is from equilibrium (see Fig. |l^). Further, Sg (< 2S^) is only slightly less than its upper 
limit 2SjP and remains so even far from equilibrium. So unlike Sg ^< LSg^, Sb is rather insensitive to the non- 
equilibrium nature of the system. Since this is coarse-grained, it is not too surprising that this is so. It would be more 
revealing to consider the behavior of fg as the size of the bins decreases, but this is currently too computationally 
intensive. 

Thermodynamic quantities, such as the entropy, allow us to investigate the underlying dynamics of the theory and 
also probe the possible deviations of the physical observables from their equilibrium values under non-equilibrium 
conditions. We shall analyze these questions below. Let us first obtain the specific heat, Cy, from the equilibrium 
ensembles; Cy may be obtained using the standard formula, 



(E 2 )eq - (E) 



Cv= - -\ 2 '- ,EQ , (49) 

where E is the energy per site. We find that the specific heat has a weak temperature dependence which may be fitted 
by a simple temperature dependence, 

CV = C T- Q , C* = 0.86(2), a = 0.025(6) . (50) 

This entails that the energy of the lattice per site has the following behavior 



E=^-T 1 - a E = C T . i.Mi 



Since the temperature dependence is weak, it is natural to compare Cy against (E) /T obtained locally both in 
equilibrium and in non- equilibrium. Such a comparison of Cy, (E)/T for equilibrium, near equilibrium and far from 
equilibrium is shown in Fig. [l5|. We find that they all agree quite well. There seems to be an intriguing tendency 
for the non-equilibrium values of (E)/T to be larger than their equilibrium counterparts. However, this cannot be 
delineated within the errors and more investigation is necessary to see if this is a real effect. 

Let us now move onto the computation of the Gibbs entropy, Sg, in equilibrium. From the application of the first 
law of thermodynamics, we obtain 

Sg = f^dT. (52) 



T 

S G (T) = S G ,o- — (T- a -l) Q -^ 5 G , + C lnr . (53) 

In Fig. ^ we plot the computed Boltzmann entropy Sg^ as a function of temperature for systems in equilibrium, as 

well as for systems near and far from equilibrium and find a similar form, Sg = Sb,o + Sb,i logT. While we cannot 

measure Sg away from equilibrium, Sb,i — Co so that the behavior of Sg^ in equilibrium and in non- equilibrium is 

similar to that of Sg in equilibrium. We also find that Sg^ is completely insensitive to the departure from equilibrium. 
Consequently, Boltzmann's entropy and some of the thermodynamic concepts such as temperature, still retain their 
relationships in these local states far from equilibrium, even when the momentum distributions are no longer gaussian. 
Other physical quantities were also studied for effects of the system being in non-equilibrium, including T 11 , 

T 11 = \^ + \ - £ . (54) 

We find that the results for T 11 in equilibrium, near equilibrium and far from equilibrium are consistent with each 
other in a manner similar to T 00 in Fig. [L5l 
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D. Boundary Temperature Jumps 



Boundary temperature jumps in non-equilibrium steady states, such as systems experiencing thermal gradients are 
well known [(jlpfl, although it seems that their behavior has never had a suitable explanation. Systems experiencing 
shearing due to moving walls also display jumps in the velocity profile, with the fluid inside the wall moving slightly 
slower than the wall velocity, for large velocities. Such effects are known to be sources of error in experiments that 
measure transport coefficients p^j . We find that it is possible to achieve a quantitative understanding of these effects 
using simple kinetic arguments. 

We have found in this system that c s and CV are of order unity, and have at most a weak temperature dependence 
(see Fig. ^ and Fig. |l5|). The thermal conductivity is related to the mean free path by k ~ Cyc s l by a standard 
kinetic theory argument. Hence we expect k ~ £. Strictly speaking, £ is the mean free path in equilibrium, but we 
will refer to this loosely as the mean free path also away from equilibrium, to use as a natural length scale in the 
system. It corresponds to the mean free path when the thermal gradients are not too strong (up to LE-E regime in 
Table 1), but the standard kinetic theory arguments presumably break down when local equilibrium no longer holds. 
The boundary temperature jumps are due to the mean free path being non-zero and the deviation of the system from 
equilibrium. When l<L, 



on 



(55) 

boundary 



where n denotes the normal to the boundary |3jJ. This formula should apply when the jump \Tf — 7i|/Tj is relatively 
small. On dimensional grounds, the coefficient n should be on order of the mean free path £, so that 77 ~ n. We 
have verified this relation by plotting n measured directly from thermal profiles, as a function of the extrapolated 
temperature T at the boundary These are summarized in Fig. [l8] (top). In the figure we display only the jumps on 
the low temperature end since the errors are much smaller. (The large gradient data on the high temperature end, 
while consistent with the low end data, is quite noisy.) We also show a fit to the data (dashes), which gives 

ij(T) = (6.1±0.5)r- 1 - 5±0 - 1 . (56) 

On the same figure we also show the behavior of k(T), which indicates that n ~ k is consistent with the observed 
behavior. 

An independent verification of the behavior of the jumps can be made by studying how the jumps depend on the 
heat flux. We let 77 = o>k, where a is a constant to be determined. From Eq. (|55|), we can then associate the heat 
flux with the right side of Eq. ( |55| ) : 

Ti — Xf ~ a(T 01 ) NE ~ «(T 2 ° - T°)^^ + • • • . (57) 
By plotting the boundary jumps directly as a function of the heat flux in non-equilibrium steady states both near 



and far from equilibrium, we obtain the data in Fig. 18 (bottom). One can see that there is a simple relationship, 
where the slope gives a — 2.6(1), consistent with the assumption that n = an. The understanding of these jumps 
together with that of the temperature profile (|39|) provides a complete description of T(x) in terms of the boundary 



temperatures, (T-f^T-P). As was alluded to earlier, these simple relations (55),(p7|) break down in the LNE region of 
Table 1, where the thermal gradients are too strong. 

A note on boundary conditions is in order: There is reasonable freedom in how the thermostats are implemented. 
We can vary the number of thermostatted sites or how the demons are couple to the physical degrees of freedom. We 
find that the changes in the way thermostats are implemented bring about changes in the boundary jumps when we 
are sufficiently far from equilibrium. This will also result in a different (T 01 ) in a manner consistent with Eq. J4C|). 
Different thermostats will correspond to different values of a, which is not an intrinsic parameter of the </> 4 theory 
but rather reflects different manners in which a heat bath might efficiently couple to the scalar field theory. In other 
words, as expected, different thermostats do not change our understanding of the physics at all. 

In understanding the generality of Eq. ( |57j ) , we note that the FPU (3 model ( p8[ ) displays temperature profiles that 
depend on L. In Refs. (g), it has been shown that the heat flux varies with the size as (T 01 ) ~ L~ 55 . Consequently, 
one would expect that the boundary jumps seen there behave as 5T ~ L~ 055 . A cursory analysis of those thermal 
profiles suggests that this is indeed the case. 
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E. Non-Equilibrium Distributions 



It is worth making a few remarks on the non-equilibrium distribution functions we have seen. Many usual approaches 
to the non-equilibrium statistical properties of systems make from the outset certain model dependent choices. In 
some thermo-field approaches, the temperature profile T(x) is chosen to have a specific form, resulting in a statistical 
operator of the type exp(— H(ir, <fi)/T(x)), where H(ir, 4>) might include contributions due to transport mM. Employing 
the approach of Zubarev Jl6{ , one assumes a particular behavior for the steady state distribution, which is then used 
to solve the dynamical equations. 

The local non-equilibrium distribution, fg\(j>k>'^k) is shown in Fig. 17 for a system very far from equilibrium. 



While the distribution / in the full configuration space is expected to be fractal, fg\(t>k,^k) is only a projection 
on a lower dimensional space, and should be smooth. Such a non-equilibrium statistical operator will differ from 
those assumed in non-equilibrium approaches in their correlations. One aspect of the fractal nature of / is that 
the dynamical space is of reduced dimension, so that additional correlations will exist that are not present in those 
standard non-equilibrium approaches. 

Another limit in which one obtains unusual thermal profiles is a 'ballistic' limit in which the mean free path is on 
the order or larger than the size of the system. Because the excitations pass through the system rather readily, one 
has large boundary jumps on both ends. As a consequence, the thermal profile is almost flat, in spite of the values of 
the endpoint temperatures. In this case the system is still thermalized; the distributions of momenta at the endpoints 
and inside are gaussian. These types of results are included in our figures of entropy, specific heat and so forth, in 
the near-equilibrium data set, and they fall in line with non-ballistic results. 



F. Local Equilibrium 

While we do not yet have any precise criteria regarding when local equilibrium fails to be a good approximation, we 
can analyze various physical quantities and compare them to their values when local equilibrium holds. This should 
provide us with a measure of the dependence of the physics on local equilibrium. While field dependent observables, 
such as moments of <j>(x) will depend on the model under investigation, the momenta should be more robust. A 
natural place to start is with the analysis of cumulants. If local equilibrium is expected holds, we should have 

(4)=T, <^)=3T 2 , (4) = 1ST 3 (58) 

and so forth. Equivalently, we can say that 

{{<)) = (4) - 3(^) 2 , ((4)) = (4) lb{4)(4) + 30(^} 3 . (59) 

vanish in local equilibrium. Hence one measure of the breaking of local equilibrium when we are far from equilibrium 
is 



3(n 2 ) 2 3(7^ 



(60) 



Similar expressions for higher cumulants are also possible. 

Non-equilibrium behavior begins to emerge in certain quantities as soon as T\ ^ T<i. For instance, 
(tt (x)} / (tt 2 (x)) 2 > 3, the equality holding in equilibrium, and the value growing as one departs from equilibrium. One 
can also see that the non-equilibrium measure is not locally Boltzmann since quantities such as (tt(x)4>(x')} ^ for 
x ^ x\ as is expected for a system supporting some type of transport. However, the origins in this case we attribute 
not to additional terms in the statistical distribution /, but rather to additional correlations in non-equilibrium mea- 
sure due to its reduced dimensionality, as expected from dynamical systems theory. Unfortunately, there is no way 
to explicitly estimate the dimensional loss without direct measurement of the entire Lyapunov spectrum, which is a 
rather numerically intensive task. 

As one moves away from equilibrium, the linear response theory is expected to eventually breakdown. Such behavior 
is confirmed for our system; when the temperature gradient becomes too large, the formula for the heat flow ([it)]) ceases 
to be valid. We display the relative difference of the measured current to the current obtained from the linear response 
theory ( f40| ) in Fig. |l^ (bottom). It can be seen that the relative deviations can be of order one for large thermal 
gradients signaling the breakdown of the linear response theory. The deviations are plotted against k(T)VT/T which 
seems to be the natural scale, since k is the mean free path, roughly speaking, which is the natural length scale in 
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the problem as discussed in section |y| D. One obvious possibility would be to interpret this as a nonlinear response of 
the thermal conductivity, k(T, T 01 ), which may, for instance, be parametrized as 

«(T, T 01 ) = «o(T) + k 2 (T) (T 01 ) 2 + O ((T 01 ) 4 ) . (61) 



Such approaches have been discussed in the literature 17 



However, such an interpretation presupposes, perhaps tacitly, that the concept of local equilibrium holds and that 
the standard notion of temperature applies, amongst other things. Therefore, it is imperative to first check that the 
local equilibrium is achieved in this 'non-linear' regime and that it is this point that we now wish to investigate with 
care. Such questions have been asked previously and in those situations, the local equilibrium was seen to be valid 
p5[. First, we need to ask ourselves what constitutes local equilibrium? The concept of local equilibrium has been 



defined previously [[16 36 and it is not our intention here to discern the possible differences in the various definitions 
of local equilibrium. In the least, local equilibrium assumes that we have a Maxwellian distribution for the momentum 
for the class of Hamiltonians we work with, leading to the usual concept of temperature, which will be our point of 
investigation. 

In Fig. |l^ (top), we plot the fourth cumulant of 7T, (|60|), against k(T)VT/T. The cumulant, which is defined 
locally, quite clearly deviates from the equilibrium value under strong thermal gradients. Even in these situations, 
the thermostatted boundary sites are in local equilibrium, as they should be, as shown in Fig. It can be seen 
that the deviations from local equilibrium start to occur at roughly the same value of k(T)VT/T ~ 1/10 where linear 
response theory breaks down. In Table 1, we identify this as a steady state which is locally non-equilibrium (LNE). 
We have verified that higher cumulants display similar behavior. Therefore, at least in our model, the breakdown of 
local equilibrium needs to be considered when nonlincarity of the response is to be analyzed. Of course, this does not 
preclude the possibility of the nonlinearity of the response as discussed above, but that the nonlinearity needs to be 
disentangled from the deviations from local equilibrium with care. 

We emphasize here that a priori, this needs not be the case. Namely, it is in principle possible that there is a 
region where the concept of local equilibrium is still valid and yet the linear response theory breaks down. In such a 
situation, 'non-linear response' theory should be quite appropriate for analyzing the situation. However, in the cases 
we studied, such regimes do not exist and the failure of the linear response theory occurs simultaneously with the 
breakdown of local equilibrium. 



VI. CONCLUSIONS 



We have constructed non-equilibrium steady states for classical 4 lattice field theory in one dimension, under 
conditions near and far from equilibrium. We obtained the behavior of the thermal conductivity with respect to 
the temperature in the linear regime and found that the results were consistent with the linear response regime. 
The underlying dynamics of the theory was investigated and physical quantities such as the speed of sound, heat 
capacity, Boltzmann's entropy and their temperature dependence, were obtained. The results could consistently 
understood using the kinetic theory approach. This understanding was further used to clarify the dynamics behind 
the temperature gaps that arise at the boundaries of the system. We also found that for temperature gradients that 
are not too large, the linear response law is adequate for understanding the behavior of the system, even though the 
temperature profile might be visibly non-linear. For even larger gradients, even though the system is in a steady state, 
we found that the linear response law eventually ceases to hold, but the local equilibrium is also violated. 

We have classified the steady states in Table 1, which identify distinct dynamical regimes of the theory. It would be 
nice to develop more precise measures for these dynamical regimes, but it is clear that even a one component classical, 
lattice field theory does contain a means to understand the non-equilibrium physics of many-body systems. It would 
be interesting to extend these results to theory with phase transitions, multi-component theories which would also 
allow the analysis of Onsager reciprocity relations and so on. The additional degrees of freedom will provide additional 
measurable quantities, but we expect many of the qualitative features of this simple model to persist. 
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FIG. 1. Comparison between equilibrium (left) and non-equilibrium (right) steady state distribution functions for momenta 
-Kk and heat flux T 01 . The measured temperature at the site k is indicated by T, and the boundary temperatures by Tj , T° . All 
measurements are made at mid lattice, x = k = L/2, with the exception of the lower right figure which was made at x — L/4. 
The heat flux (T 01 ) is zero in equilibrium, so f(T 01 ) is symmetric. Away from equilibrium, (T 01 ) 7^ 0, so the distribution 
develops asymmetry as seen in the lower right. 
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FIG. 2. Time evolution of the heat flux (T 01 )ne at three points inside a system with L = 800, denoted by the ratio x/L, 
and (J?, T°) = (0.3,0.7). As there are no sources or sinks inside, the asymptotic values must converge to the same value and 
they do. 
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FIG. 3. Finite temperature autocorrelation functions for T^ v as a function of time scaled by the mean free time r ~ k(T). 
We plot the absolute value of C(t) = ((T"" (t)T^ (0)) - {T M ' y (0)> 2 )/{(T' ll '(0)) 2 ). At large times, these functions oscillate about 
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FIG. 4. Green-Kubo integrals n(T,t) up to time t for lattices with L — 100. As t — > oo, n(T,t) — > k(T). We plot the ratio 
of n(T,t) to its asymptotic value versus time, normalized by the mean free time r. Since Cy ~ c s ~ 1, r is approximated by 
the magnitude of the conductivity at that temperature. The dashed lines are the anticipated behaviors if the long-time tail 
divergences were present. One can see that on the time scales up to t ~ lOr, such behavior can be seen, although it vanishes 
for larger times. 
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FIG. 5. Power law behavior of the thermal conductivity n for lattice <f> i theory obtained from near equilibrium (□) and 
Green-Kubo (x) measurements for various lattice sizes L. The power law fit k(T) = 2.83(4)/T 1 ' 35(2) is indicated by the dashes. 
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FIG. 6. Spatial dependence of the equilibrium auto-correlation function G(x, t; 0, 0) at selected times, x = is the center of 
the lattice. By examining the peaks, one can estimate how fast the excitations propagate through the system. Here L — 160 
and T = 1/10. 
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FIG. 7. Temperature dependence of the "speed of sound", c s . The measurements were made in a system with L = 163 but 
the speed was also seen to be independent of the size. 
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FIG. 8. A typical near equilibrium temperature profile (solid), for L = 8000, compared to a linear profile (dashes). The two 
endpoint thermostats are TJ 3 = 0.15 and T° = 0.35. 
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FIG. 9. Method used for direct computation of the thermal conductivity near equilibrium. In this regime the thermal profiles 
are linear and Fourier's law is well reproduced. By taking increasingly small differences T2 — if around a fixed temperature 
T, one can measure (T 01 )ne for many values of VT and extract the slope, which equals the conductivity. Here we show such 
a result for L = 800 at T = 1/4. 
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FIG. 10. Evidence of bulk behavior of the thermal conductivity for selected temperatures T = 1 (O), T = 1/2 (x), T = 1/4 
(□) and T = 1/10 (+). The dashed lines are the predictions from the power law fit (|3q) to k(T). 
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FIG. 11. As the boundary temperature differences increase, the system departs from the constant gradient regime and the 
thermal profiles develop curvature. In the figure we show, the boundary temperatures are (T®, ) = (0.3, 0.7) (dots), (0.2, 0.8) 
(dashes) and (0.1,2.0) (solid) for a system with L — 800. 
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FIG. 12. (Top) Non-equilibrium thermal profile for boundary temperatures {T^,T^) = (0.05,1) (solid) compared to pre- 
dictions (solid). One can see deviations on the of a few %, most notably at the high temperature end (bottom, right). The 
boundary jumps are about equal and are readily understood. 
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FIG. 13. As the system is driven even further from equilibrium, we begin to see departures from local equilibrium and linear 
response. Here we show the thermal profile for a system with T$ /T° = 100 and L — 163 (solid), compared to the theory curve 
from Eq. (^) (dashes). The agreement is not good. 
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FIG. 14. (Top) Steady state distribution of momenta itk in the center of the strongly non- equilibrium system discussed in 
the previous figure. We observe that these distributions are characteristically more strongly peaked than the local equilibrium 
gaussian at that temperature (solid). (Bottom) The ratio of the theory profile to the measured profile in the previous figure, 
as a function of position. The agreement is worse for lower temperatures. 
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FIG. 15. Specific heat, Cv, measured in equilibrium is plotted as a function of temperature. The ratio of the energy density 
(T 00 ) is also plotted against the temperature T(x), at equilibrium (O), near equilibrium (+) and far from equilibrium (□). We 
see that there is a temperature dependence and as well as seemingly a slight effect due to non-equilibrium physics. 
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FIG. 16. Coarse-grained one-body Boltzmann entropy Sb(T) as a function of T for equilibrium (O), near equilibrium (+) 
and far from equilibrium (□) systems. One can see that the (local) entropy is insensitive to the departure from equilibrium. 
The solid line indicates a logarithmic fit, Sb(T) = Sb,o + Sb,i logT, with Sb,o = 2.5(1) and Sb,i = 0.80(2). 
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FIG. 17. (Top) Non-equilibrium steady-state distribution function f(n k ,(f> k ) for a system with {T®,T$) = (0.1,6), where 
k{= L/2) is taken to be the middle of the system. 
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FIG. 18. (Top) Behavior of boundary jumps as a function of temperature. r/(T) can be seen to obey a power law behavior 
(dashes) similar to the conductivity. (Bottom) Evidence that the boundary jumps are simply related to the heat flux. 
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FIG. 19. Measures of the breaking of local equilibrium. (Top) The behavior of the cumulants of the momentum distribution 
for various non-equilibrium steady states, as a function of kVT/T ~ £VT/T, where £ is the mean free path. One can see 
that noticeable departures from gaussian momentum distributions develop for nVT/T > 1/10. (Bottom) Departure of the 
predicted heat flux from the measured heat flux. One can see that the predictions, based on local equilibrium and the shape 
of the thermal profile, begin to deviate from measured values in the same regime. 
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FIG. 20. Behavior of the fourth moment (% 4 ) normalized by (vr 2 ) 2 . (Top) We show the time evolution at the center of a 
system with L = 160 sites and boundary temperatures (Tf,T%) — (0.1, 1). We see that the ratio is well converged. (Bottom) 
We show the behavior of the moments along the lattice. The local equilibrium value is indicated by the dashed line. While 
the endpoints are in thermal equilibrium one can see that the low temperature end suffers the largest departures from local 
equilibrium. 
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